%% Aim
% Conduct PTM modeling

%% Setup
clc;
clear;
tic;
root_path = pwd; 
addpath(genpath(root_path));

nrole=2;
nprice=10;
nlottery=10;

role_label={'Buyer','Seller'};
role_color_dark={[235/255 112/255 99/255],[106/255 142/255 208/255]}; 
role_color_median={[238/255 137/255 126/255],[144/255 171/255 220/255]}; 
role_color_light={[243/255 171/255 163/255],[165/255 187/255 227/255]};

bias_label={'Val','Res','Val vs. Res'};
bias_color_median={[195/255 205/255 147/255],[192/255 173/255 203/255]}; 

%% Load data
beh_data=readtable([root_path,'\data\behavior\behavior.csv']);
sub_data_raw=readtable([root_path,'\data\subject\subject.csv']);
sub_data_beh=readtable([root_path,'\data\subject\subject_behavior.csv']);
sub_data_ddm=readtable([root_path,'\data\subject\subject_ddm.csv']);
sub_data_ddm_pred=readtable([root_path,'\data\subject\subject_ddm_pred.csv']);
sub_data_gaze=readtable([root_path,'\data\subject\subject_gaze.csv']);
sub_data_pupil=readtable([root_path,'\data\subject\subject_pupil.csv']);

%% Incorporate data
sub_data=[sub_data_raw,sub_data_pupil(:,4:end),sub_data_ddm_pred(:,end-1:end)];
nsub=height(sub_data);

%% 
        









%% PTM1

% params(1): lottery weighting for buyer
% params(2): decisiveness
% params(3): response bias for buyer
% params(4): lottery weighting for seller
% params(5): response bias for seller
% params(6): fixed utility for buyer
% params(7): fixed utility for seller

idx=1;
for isub = 1:height(sub_data)

    theta = [1    1    0    1    0    0    0 ]; % Initial values
    lb =    [0    0   -1    0   -1  -Inf -Inf];
    ub =    [Inf Inf   1   Inf   1   Inf  Inf];    
    options = optimset('Display','iter','TolFun',1e-8,'MaxFunEvals',10000,'Display','notify','MaxIter',10000);   
    wdata = beh_data(beh_data.Subject==sub_data.Subject(isub),:);
    [MLE{idx}(isub,1:length(theta)),NegLogL(isub,idx),exitflag(isub,idx)] = fmincon(@ptm1,theta,[],[],[],[],lb,ub,[],options,wdata);
    LogL(isub,idx) = -NegLogL(isub,idx);
    [AIC(isub,idx),BIC(isub,idx)] = aicbic(LogL(isub,idx),length(theta),height(wdata));
    
end

sub_data.PTMwLotteryBuy=MLE{idx}(:,1);
sub_data.PTMwLotterySell=MLE{idx}(:,4);

sub_data.PTMzBuy=MLE{idx}(:,3);
sub_data.PTMzSell=MLE{idx}(:,5);

sub_data.PTMInterBuy=MLE{idx}(:,6);
sub_data.PTMInterSell=MLE{idx}(:,7);

sub_data.PTMMiuBuy=MLE{idx}(:,2);
sub_data.PTMMiuSell=MLE{idx}(:,2);

sub_data.PTMResBias=sub_data.PTMzSell-sub_data.PTMzBuy;
sub_data.PTMValBias=sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy;
sub_data.PTMInterBias=sub_data.PTMInterSell-sub_data.PTMInterBuy;
sub_data.PTMMiuBias=sub_data.PTMMiuSell./sub_data.PTMMiuBuy;

sub_data.PTMValBiasLog=log(sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy);

% table
A=mean([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
B=std([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
rowNames = {'mean','std'};
colNames = {'PTMwLotteryBuy','PTMwLotterySell','PTMValBias','PTMzBuy','PTMzSell','PTMResBias','PTMInterBuy','PTMInterSell','PTMInterBias','PTMMiuBuy','PTMMiuSell','PTMMiuBias','PTMValBiasLog'};
disp(array2table([A;B],'RowNames',rowNames,'VariableNames',colNames));
[A;B];

%% PTM1 predict

sdata=sub_data;

temp =table;
temptemp=table;
for isub=1:height(sdata)
    bdata = beh_data(beh_data.Subject==sdata.Subject(isub),:);
    temptemp = (bdata.Role==1).*(...
                                 (sdata.PTMzBuy(isub)>=0 ).*( (1-sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))+sdata.PTMzBuy(isub) ) + ...
                                 (sdata.PTMzBuy(isub)< 0 ).*  (1+sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))...
                                 ) + ...
               (bdata.Role==2).*(...
                                 (sdata.PTMzSell(isub)>=0).*( (1-sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))+sdata.PTMzSell(isub)) + ...
                                 (sdata.PTMzSell(isub)< 0).*  (1+sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))...
                                 ) ;
    temp=[temp;array2table(temptemp)];
end

beh_data=beh_data;
beh_data.PTMProbPred=temp.temptemp;

beh_data.PTMEndingPred=beh_data.PTMProbPred>=0.5;
beh_data.PTMChoicePred = (2-beh_data.Role) .* beh_data.PTMEndingPred + (beh_data.Role-1) .* (1-beh_data.PTMEndingPred);

for isub=1:height(sub_data)
    bdata=beh_data(beh_data.Subject==isub,:);
    sub_data.PTM1ChoicePred(isub)=mean(bdata.Choice==bdata.PTMChoicePred); 
    ChoicePred(isub,idx)=mean(bdata.Choice==bdata.PTMChoicePred); 
end   

%% Prediction Accuracy compared with DDM1

mean(sub_data.PTM1ChoicePred)
std(sub_data.PTM1ChoicePred)
mean(sub_data.DDMChoicePred)
std(sub_data.DDMChoicePred)

mean(sub_data.PTM1ChoicePred-sub_data.DDMChoicePred);
std(sub_data.PTM1ChoicePred-sub_data.DDMChoicePred);

[h,p,ci,stats]=ttest(sub_data.PTM1ChoicePred-sub_data.DDMChoicePred, 0);%,'tail','right');
[h,p,ci,stats]=ttest(sub_data.DDMChoicePred, sub_data.PTM1ChoicePred);%,'tail','right');

% histogram
figure('Renderer', 'painters', 'Position', [10 10 370 350]);
ax=subplot(1,1,1);
%histogram(sub_data.PTM1ChoicePred,0.7:0.02:1,'EdgeColor',[0.5 0.5 0.5],'FaceColor',[0.5 0.5 0.5]);
histogram(sub_data.PTM1ChoicePred,0.7:0.02:1,'EdgeColor',role_color_median{2},'FaceColor',role_color_median{2});
hold on
%histogram(sub_data.DDMChoicePred,0.7:0.02:1,'EdgeColor',[0 0 0],'FaceColor',[0 0 0]);
histogram(sub_data.DDMChoicePred,0.7:0.02:1,'EdgeColor',role_color_median{1},'FaceColor',role_color_median{1});
set(gca,'TickDir','out');
xlabel('Accuracy');
ylabel('Number of Participants');
%xlim([0.65 1]);
xticks(0.7:0.1:1.0);
%yticks(0:0.1:0.4);
print(1, '-dtiff', [root_path, '\figure\model\ptm\ptm_ddm_accuracy.tif'], '-r200');
hold off
close;

%% PTM1 predict: Plot

    for irole=1:nrole
        for iprice=1:nprice
            for ilottery=1:nlottery
                bdata=beh_data(  beh_data.Role==irole & ...
                                 beh_data.Price==iprice & ...
                                 beh_data.Lottery==ilottery*2,:);
                prob_accept.matrix.group.role(irole).value(ilottery,iprice)=nanmean(bdata.PTMChoicePred);
            end
        end
    end

    figure('Renderer', 'painters', 'Position', [10 10 1000 350]);
    for irole=1:nrole
        subplot(1,2,irole);
        A=prob_accept.matrix.group.role(irole).value;
        h=imagesc(A);
        cmap_up=[linspace(1,192/256,1000)',linspace(1,0,1000)',linspace(1,0,1000)']; 
        cmap_down=[linspace(0/256,1,1000)',linspace(32/256,1,1000)',linspace(102/256,1,1000)']; 
        cmap=[cmap_down;cmap_up];
        colormap(cmap);
        set(gca,'TickDir','out');
        xlabel('Price');
        ylabel('Lottery');
        yticks(1:nprice);
        xticks(1:nlottery);
        yticklabels({'2','4','6','8','10','12','14','16','18','20'});
        xticklabels({'1','2','3','4','5','6','7','8','9','10'});
        colorbar;
        caxis([0 1]);
        title(role_label{irole});
        hold on
    end
    hold off
    print(1, '-dtiff', [root_path, '\figure\behavior\choice\matrix\prob_accept_matrix_group_pred_ptm.tif'], '-r200');
    close;       

%%










%% PTM1: Plot subject estimates

sdata=sub_data;

figure('Renderer', 'painters', 'Position', [10 10 380 350]);
ax=subplot(1,1,1);
x1=sdata.PTMwLotteryBuy;
y1=sdata.PTMzBuy;
[r1,p1]=corr(x1,y1);
scatter(x1,y1,36,role_color_median{1},'filled');
hold on
x2=sdata.PTMwLotterySell;
y2=sdata.PTMzSell;
[r2,p2]=corr(x2,y2);
scatter(x2,y2,36,role_color_median{2},'filled');
ylim([-0.3 0.3]);
yticks(-0.3:0.1:0.3);
xlim([0 1.2]);
set(ax,'TickDir','out');
ylabel('Resopnse tendency');
xlabel('Weighting of Lottery');
box on;
line(xlim,[0 0],'Color',[0 0 0],'LineStyle','--');
print(1, '-dtiff', [root_path,'\figure\model\ptm\ptm_estimates.tif'], '-r200');
hold off;
close;

%% Plot PTM biases log

sdata=sub_data;

figure('Renderer', 'painters', 'Position', [10 10 380 350]);
ax=subplot(1,1,1);
x=sdata.PTMValBiasLog;
y=sdata.PTMResBias;
[r,p]=corr(x,y);
scatter(x,y,36,[0.75,0.75,0.75],'filled');
hold on
set(ax,'TickDir','out');
ylabel('PTM Response Bias');
xlabel('PTM Valuation Bias');
box on;

%ylim([-0.2 0.3]);
%xlim([0 3]);

ylim([-0.1 0.2]);
yticks(-0.1:0.1:0.2);
line(xlim,[0 0],'Color',bias_color_median{2},'LineStyle','--');
line([0 0],ylim,'Color',bias_color_median{1},'LineStyle','--');

print(1, '-dtiff', [root_path,'\figure\model\ptm\ptm_bias_log.tif'], '-r200');
hold off;
close;

%% PTMRes, PTMVal and PTM Intercept: T Test

% PTM Val
sdata=sub_data;
mean(sdata.PTMValBiasLog)
std(sdata.PTMValBiasLog)
[h,p,ci,stats]=ttest(sdata.PTMValBiasLog,0);%,'tail','right');

% PTM Res
mean(sdata.PTMResBias)
std(sdata.PTMResBias)
[h,p,ci,stats]=ttest(sdata.PTMResBias,0);%,'tail','right');
sdata=sub_data;
mean(sdata.PTMzBuy)
std(sdata.PTMzBuy)
[h,p,ci,stats]=ttest(sdata.PTMzBuy,0);%,'tail','right');
mean(sdata.PTMzSell)
std(sdata.PTMzSell)
[h,p,ci,stats]=ttest(sdata.PTMzSell,0);%,'tail','right');

% PTM Intercept
sdata=sub_data;
mean(sdata.PTMInterBuy)
std(sdata.PTMInterBuy)
[h,p,ci,stats]=ttest(sdata.PTMInterBuy,0);%,'tail','right');
mean(sdata.PTMInterSell)
std(sdata.PTMInterSell)
[h,p,ci,stats]=ttest(sdata.PTMInterSell,0);%,'tail','right');
mean(sdata.PTMInterBias)
std(sdata.PTMInterBias)
[h,p,ci,stats]=ttest(sdata.PTMInterBias,0);%,'tail','right');

%% Quadrant

% DDM
mean(sub_data.ValBiasLog>0 & sub_data.ResBias>0);
mean(sub_data.ValBiasLog<0 & sub_data.ResBias>0);
mean(sub_data.ValBiasLog>0 & sub_data.ResBias<0);
mean(sub_data.ValBiasLog<0 & sub_data.ResBias<0);

% Similarity between DDM and PTM
q1=height(sub_data(sub_data.PTMValBiasLog>0 & sub_data.ValBiasLog>0 & sub_data.PTMResBias>0 & sub_data.ResBias>0,:));
q2=height(sub_data(sub_data.PTMValBiasLog<0 & sub_data.ValBiasLog<0 & sub_data.PTMResBias>0 & sub_data.ResBias>0,:));
q3=height(sub_data(sub_data.PTMValBiasLog<0 & sub_data.ValBiasLog<0 & sub_data.PTMResBias<0 & sub_data.ResBias<0,:));
q4=height(sub_data(sub_data.PTMValBiasLog>0 & sub_data.ValBiasLog>0 & sub_data.PTMResBias<0 & sub_data.ResBias<0,:));
% Overall
(q1+q2+q3+q4)/64;
% Similarity by bias
height(sub_data((sub_data.PTMValBiasLog>0 & sub_data.ValBiasLog>0) | (sub_data.PTMValBiasLog<0 & sub_data.ValBiasLog<0),:))/64;
height(sub_data((sub_data.PTMResBias>0 & sub_data.ResBias>0) | (sub_data.PTMResBias<0 & sub_data.ResBias<0),:))/64;

% for reviewer
q1/sum(sub_data.ValBiasLog>0 & sub_data.ResBias>0);
q2/sum(sub_data.ValBiasLog<0 & sub_data.ResBias>0);
q3/sum(sub_data.ValBiasLog<0 & sub_data.ResBias<0);
q4/sum(sub_data.ValBiasLog>0 & sub_data.ResBias<0);


%%










%% Correlation between PTM and DDM Biases

sdata=sub_data;

[r,p]=corr([sdata.PTMValBiasLog,sdata.PTMResBias,sdata.PTMInterBias,sdata.ValBiasLog,sdata.ResBias,sdata.vInterceptBias,sdata.aBias]);

    figure('Renderer', 'painters', 'Position', [10 10 450 360]);
        subplot(1,1,1);
        A=r;
        h=imagesc(A);
        cmap_up=[linspace(1,192/256,1000)',linspace(1,0,1000)',linspace(1,0,1000)']; 
        cmap_down=[linspace(0/256,1,1000)',linspace(32/256,1,1000)',linspace(102/256,1,1000)']; 
        cmap=[cmap_down;cmap_up];
        colormap(cmap);
        set(gca,'TickDir','out');
        %xlabel('Price');
        %ylabel('Lottery');
        yticks(1:nprice);
        xticks(1:nlottery);
        yticklabels({'2','4','6','8','10','12','14','16','18','20'});
        xticklabels({'1','2','3','4','5','6','7','8','9','10'});
        colorbar;
        caxis([-1 1]);
        hold on
    hold off
    print(1, '-dtiff', [root_path, '\figure\model\ptm\ptm_ddm_post_corr.tif.tif'], '-r200');
    close;     


%% PTM biases vs. gaze

sdata=sub_data;

% Normalize data for beta size comparison
sdata.PTMValBias_scaled=rescale(sdata.PTMValBias);
sdata.PTMResBias_scaled=rescale(sdata.PTMResBias);
sdata.PTMValBiasLog_scaled=rescale(sdata.PTMValBiasLog);
sdata.Role1_scaled=rescale(sdata.Role1);

% Main
glm = fitglm(sdata,'GazeDurLotteryRatioBiasMinus ~ 1 + PTMValBiasLog_scaled ');
glm = fitglm(sdata,'GazeDurLotteryRatioBiasMinus ~ 1 + PTMResBias_scaled ');
glm = fitglm(sdata,'GazeDurLotteryRatioBiasMinus ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled ');
glm = fitglm(sdata,'GazeDurLotteryRatioBiasMinus ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled + Role1');

% Control: val bias as ratio
glm = fitglm(sdata,'GazeDurLotteryRatioBiasMinus ~ 1 + PTMValBias_scaled ');
glm = fitglm(sdata,'GazeDurLotteryRatioBiasMinus ~ 1 + PTMResBias_scaled ');
glm = fitglm(sdata,'GazeDurLotteryRatioBiasMinus ~ 1 + PTMValBias_scaled + PTMResBias_scaled ');
glm = fitglm(sdata,'GazeDurLotteryRatioBiasMinus ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled + Role1');

% Control: gaze bias as ratio
glm = fitglm(sdata,'GazeDurLotteryRatioBias ~ 1 + PTMValBiasLog_scaled ');
glm = fitglm(sdata,'GazeDurLotteryRatioBias ~ 1 + PTMResBias_scaled ');
glm = fitglm(sdata,'GazeDurLotteryRatioBias ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled ');
glm = fitglm(sdata,'GazeDurLotteryRatioBias ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled + Role1');

% Control: gaze bias as ratio and val bias as ratio
glm = fitglm(sdata,'GazeDurLotteryRatioBias ~ 1 + PTMValBias_scaled ');
glm = fitglm(sdata,'GazeDurLotteryRatioBias ~ 1 + PTMResBias_scaled ');
glm = fitglm(sdata,'GazeDurLotteryRatioBias ~ 1 + PTMValBias_scaled + PTMResBias_scaled ');
glm = fitglm(sdata,'GazeDurLotteryRatioBias ~ 1 + PTMValBias_scaled + PTMResBias_scaled + Role1');

% Plot
glm = fitglm(sdata,'GazeDurLotteryRatioBiasMinus ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled ');
disp('Predicted by HDDM parameters');
disp(glm);

B=0;
SE=0;
Lower=0;
Upper=0;
% extract estimates
for i=1:3
   B(i,1)=glm.Coefficients.Estimate(i);
   SE(i,1)=glm.Coefficients.SE(i); 
end
CI=coefCI(glm);
Lower=CI(:,1);
Upper=CI(:,2);
disp('Confidence interval');
disp(CI);

% bar plot
ax=subplot(1,1,1);
set(ax, 'position', [0.15, 0.15, 0.75*1024/1280, 0.8]);
ab=bar([1 2],[B(3,1),B(2,1)],'EdgeColor','none','barwidth',0.5);
xlim([0.2 2.8]);
ab.FaceColor='flat';
ab.CData(1,:) = bias_color_median{1};
ab.CData(2,:) = bias_color_median{2};
hold on
er=errorbar([1 2],[B(3,1),B(2,1)],[SE(3,1),SE(2,1)]);
er.Color = [0 0 0];                            
er.LineStyle = 'none';  
xticklabels({'Val','Res','Inter','a'});
set(gca,'TickDir','out');
xlabel('(scaled)');
ylabel('Beta');
ylim([-0.1 0.4]);
yticks(-0.1:0.1:0.4);
%yticks(-0.4:0.4:0.8);
hold off

print(1, '-dtiff', [root_path,'\figure\model\ptm\gldr_regress_on_ptm_biases.tif'], '-r200');
close;

%% PTM biases vs. First gaze

sdata=sub_data;

% Normalize data for beta size comparison
sdata.PTMValBias_scaled=rescale(sdata.PTMValBias);
sdata.PTMResBias_scaled=rescale(sdata.PTMResBias);
sdata.PTMValBiasLog_scaled=rescale(sdata.PTMValBiasLog);
sdata.Role1_scaled=rescale(sdata.Role1);

% Main
glm = fitglm(sdata,'FirstGazeLotteryRatioBiasMinus ~ 1 + PTMValBiasLog_scaled ');
glm = fitglm(sdata,'FirstGazeLotteryRatioBiasMinus ~ 1 + PTMResBias_scaled ');
glm = fitglm(sdata,'FirstGazeLotteryRatioBiasMinus ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled ');
glm = fitglm(sdata,'FirstGazeLotteryRatioBiasMinus ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled + Role1');

% Control: val bias as ratio
glm = fitglm(sdata,'FirstGazeLotteryRatioBiasMinus ~ 1 + PTMValBias_scaled ');
glm = fitglm(sdata,'FirstGazeLotteryRatioBiasMinus ~ 1 + PTMResBias_scaled ');
glm = fitglm(sdata,'FirstGazeLotteryRatioBiasMinus ~ 1 + PTMValBias_scaled + PTMResBias_scaled ');
glm = fitglm(sdata,'FirstGazeLotteryRatioBiasMinus ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled + Role1');

% Plot
glm = fitglm(sdata,'FirstGazeLotteryRatioBiasMinus ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled ');
disp('Predicted by HDDM parameters');
disp(glm);

B=0;
SE=0;
Lower=0;
Upper=0;
% extract estimates
for i=1:3
   B(i,1)=glm.Coefficients.Estimate(i);
   SE(i,1)=glm.Coefficients.SE(i); 
end
CI=coefCI(glm);
Lower=CI(:,1);
Upper=CI(:,2);
disp('Confidence interval');
disp(CI);

% bar plot
ax=subplot(1,1,1);
set(ax, 'position', [0.15, 0.15, 0.75*1024/1280, 0.8]);
ab=bar([1 2],[B(3,1),B(2,1)],'EdgeColor','none','barwidth',0.5);
xlim([0.2 2.8]);
ab.FaceColor='flat';
ab.CData(1,:) = bias_color_median{1};
ab.CData(2,:) = bias_color_median{2};
hold on
er=errorbar([1 2],[B(3,1),B(2,1)],[SE(3,1),SE(2,1)]);
er.Color = [0 0 0];                            
er.LineStyle = 'none';  
xticklabels({'Val','Res','Inter','a'});
set(gca,'TickDir','out');
xlabel('(scaled)');
ylabel('Beta');
ylim([-0.1 0.7]);
%yticks(-0.1:0.1:0.4);
%yticks(-0.4:0.4:0.8);
hold off

print(1, '-dtiff', [root_path,'\figure\model\ptm\first_gaze_regress_on_ptm_biases.tif'], '-r200');
close;

%% PTM biases vs. pupil

sdata=sub_data;

% Normalize data for beta size comparison
sdata.PTMValBias_scaled=rescale(sdata.PTMValBias);
sdata.PTMResBias_scaled=rescale(sdata.PTMResBias);
sdata.PTMValBiasLog_scaled=rescale(sdata.PTMValBiasLog);
sdata.Role1_scaled=rescale(sdata.Role1);

% Main
glm = fitglm(sdata,'DecPupilBias ~ 1 + PTMValBiasLog_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + PTMResBias_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled + Role1');

% Control: val bias as ratio
glm = fitglm(sdata,'DecPupilBias ~ 1 + PTMValBias_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + PTMResBias_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + PTMValBias_scaled + PTMResBias_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled + Role1');

% Plot
glm = fitglm(sdata,'DecPupilBias ~ 1 + PTMValBiasLog_scaled + PTMResBias_scaled ');
disp('Predicted by HDDM parameters');
disp(glm);

B=0;
SE=0;
Lower=0;
Upper=0;
% extract estimates
for i=1:3
   B(i,1)=glm.Coefficients.Estimate(i);
   SE(i,1)=glm.Coefficients.SE(i); 
end
CI=coefCI(glm);
Lower=CI(:,1);
Upper=CI(:,2);
disp('Confidence interval');
disp(CI);

% bar plot
ax=subplot(1,1,1);
set(ax, 'position', [0.15, 0.15, 0.75*1024/1280, 0.8]);
ab=bar([1 2],[B(3,1),B(2,1)],'EdgeColor','none','barwidth',0.5);
xlim([0.2 2.8]);
ab.FaceColor='flat';
ab.CData(1,:) = bias_color_median{1};
ab.CData(2,:) = bias_color_median{2};
hold on
er=errorbar([1 2],[B(3,1),B(2,1)],[SE(3,1),SE(2,1)]);
er.Color = [0 0 0];                            
er.LineStyle = 'none';  
xticklabels({'Val','Res','Inter','a'});
set(gca,'TickDir','out');
xlabel('(scaled)');
ylabel('Beta');
ylim([-0.1 0.8]);
%yticks(-0.1:0.1:0.4);
%yticks(-0.4:0.4:0.8);
hold off

print(1, '-dtiff', [root_path,'\figure\model\ptm\pupil_bias_regress_on_ptm_biases.tif'], '-r200');
close;

%% PTM vs. DDM pupil
sdata=sub_data;

% Normalize data for beta size comparison
sdata.ResBias_scaled=rescale(sdata.ResBias);
sdata.PTMResBias_scaled=rescale(sdata.PTMResBias);
sdata.Role1_scaled=rescale(sdata.Role1);

% Main
glm = fitglm(sdata,'DecPupilBias ~ 1 + PTMResBias_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + ResBias_scaled ');
glm = fitglm(sdata,'DecPupilBias ~ 1 + PTMResBias_scaled + ResBias_scaled');
glm = fitglm(sdata,'DecPupilBias ~ 1 + PTMResBias_scaled + ResBias_scaled + Role1');

%% PTM Intercept
sdata=sub_data;
[r,p]=corr(sdata.PTMInterBias, sdata.GazeDurLotteryRatioBiasMinus);
[r,p]=corr(sdata.PTMInterBias, sdata.FirstGazeLotteryRatioBiasMinus);
[r,p]=corr(sdata.PTMInterBias, sdata.DecPupilBias);

%%










%% Model comparison

%% PTM2: no role effect on weighting

% params(1): lottery weighting for buyer
% params(2): decisiveness
% params(3): response bias for buyer
% params(1): lottery weighting for seller
% params(4): response bias for seller
% params(5): fixed utility for buyer
% params(6): fixed utility for seller

idx=2;
for isub = 1:height(sub_data)

    theta = [1    1    0        0    0    0 ]; % Initial values
    lb =    [0    0   -1       -1  -Inf -Inf];
    ub =    [Inf Inf   1        1   Inf  Inf];    
    options = optimset('Display','iter','TolFun',1e-8,'MaxFunEvals',10000,'Display','notify','MaxIter',10000);   
    wdata = beh_data(beh_data.Subject==sub_data.Subject(isub),:);
    [MLE{idx}(isub,1:length(theta)),NegLogL(isub,idx),exitflag(isub,idx)] = fmincon(@ptm2,theta,[],[],[],[],lb,ub,[],options,wdata);
    LogL(isub,idx) = -NegLogL(isub,idx);
    [AIC(isub,idx),BIC(isub,idx)] = aicbic(LogL(isub,idx),length(theta),height(wdata));
    
end

sub_data.PTMwLotteryBuy=MLE{idx}(:,1);
sub_data.PTMwLotterySell=MLE{idx}(:,1);

sub_data.PTMzBuy=MLE{idx}(:,3);
sub_data.PTMzSell=MLE{idx}(:,4);

sub_data.PTMInterBuy=MLE{idx}(:,5);
sub_data.PTMInterSell=MLE{idx}(:,6);

sub_data.PTMMiuBuy=MLE{idx}(:,2);
sub_data.PTMMiuSell=MLE{idx}(:,2);

sub_data.PTMResBias=sub_data.PTMzSell-sub_data.PTMzBuy;
sub_data.PTMValBias=sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy;
sub_data.PTMInterBias=sub_data.PTMInterSell-sub_data.PTMInterBuy;
sub_data.PTMMiuBias=sub_data.PTMMiuSell./sub_data.PTMMiuBuy;

sub_data.PTMValBiasLog=log(sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy);

% table
A=mean([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
B=std([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
rowNames = {'mean','std'};
colNames = {'PTMwLotteryBuy','PTMwLotterySell','PTMValBias','PTMzBuy','PTMzSell','PTMResBias','PTMInterBuy','PTMInterSell','PTMInterBias','PTMMiuBuy','PTMMiuSell','PTMMiuBias','PTMValBiasLog'};
disp(array2table([A;B],'RowNames',rowNames,'VariableNames',colNames));
[A;B];

%% PTM1 predict

sdata=sub_data;

temp =table;
temptemp=table;
for isub=1:height(sdata)
    bdata = beh_data(beh_data.Subject==sdata.Subject(isub),:);
    temptemp = (bdata.Role==1).*(...
                                 (sdata.PTMzBuy(isub)>=0 ).*( (1-sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))+sdata.PTMzBuy(isub) ) + ...
                                 (sdata.PTMzBuy(isub)< 0 ).*  (1+sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))...
                                 ) + ...
               (bdata.Role==2).*(...
                                 (sdata.PTMzSell(isub)>=0).*( (1-sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))+sdata.PTMzSell(isub)) + ...
                                 (sdata.PTMzSell(isub)< 0).*  (1+sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))...
                                 ) ;
    temp=[temp;array2table(temptemp)];
end

beh_data=beh_data;
beh_data.PTMProbPred=temp.temptemp;

beh_data.PTMEndingPred=beh_data.PTMProbPred>=0.5;
beh_data.PTMChoicePred = (2-beh_data.Role) .* beh_data.PTMEndingPred + (beh_data.Role-1) .* (1-beh_data.PTMEndingPred);

for isub=1:height(sub_data)
    bdata=beh_data(beh_data.Subject==isub,:);
    sub_data.PTM2ChoicePred(isub)=mean(bdata.Choice==bdata.PTMChoicePred); 
    ChoicePred(isub,idx)=mean(bdata.Choice==bdata.PTMChoicePred); 
end   

%% PTM3: no role effect on response

% params(1): lottery weighting for buyer
% params(2): decisiveness
% params(3): response bias for buyer
% params(4): lottery weighting for seller
% params(3): response bias for seller
% params(5): fixed utility for buyer
% params(6): fixed utility for seller

idx=3;
for isub = 1:height(sub_data)

    theta = [1    1    0    1        0    0 ]; % Initial values
    lb =    [0    0   -1    0      -Inf -Inf];
    ub =    [Inf Inf   1   Inf      Inf  Inf];    
    options = optimset('Display','iter','TolFun',1e-8,'MaxFunEvals',10000,'Display','notify','MaxIter',10000);   
    wdata = beh_data(beh_data.Subject==sub_data.Subject(isub),:);
    [MLE{idx}(isub,1:length(theta)),NegLogL(isub,idx),exitflag(isub,idx)] = fmincon(@ptm3,theta,[],[],[],[],lb,ub,[],options,wdata);
    LogL(isub,idx) = -NegLogL(isub,idx);
    [AIC(isub,idx),BIC(isub,idx)] = aicbic(LogL(isub,idx),length(theta),height(wdata));
    
end

sub_data.PTMwLotteryBuy=MLE{idx}(:,1);
sub_data.PTMwLotterySell=MLE{idx}(:,4);

sub_data.PTMzBuy=MLE{idx}(:,3);
sub_data.PTMzSell=MLE{idx}(:,3);

sub_data.PTMInterBuy=MLE{idx}(:,5);
sub_data.PTMInterSell=MLE{idx}(:,6);

sub_data.PTMMiuBuy=MLE{idx}(:,2);
sub_data.PTMMiuSell=MLE{idx}(:,2);

sub_data.PTMResBias=sub_data.PTMzSell-sub_data.PTMzBuy;
sub_data.PTMValBias=sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy;
sub_data.PTMInterBias=sub_data.PTMInterSell-sub_data.PTMInterBuy;
sub_data.PTMMiuBias=sub_data.PTMMiuSell./sub_data.PTMMiuBuy;

sub_data.PTMValBiasLog=log(sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy);

% table
A=mean([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
B=std([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
rowNames = {'mean','std'};
colNames = {'PTMwLotteryBuy','PTMwLotterySell','PTMValBias','PTMzBuy','PTMzSell','PTMResBias','PTMInterBuy','PTMInterSell','PTMInterBias','PTMMiuBuy','PTMMiuSell','PTMMiuBias','PTMValBiasLog'};
disp(array2table([A;B],'RowNames',rowNames,'VariableNames',colNames));
[A;B];

%% PTM1 predict

sdata=sub_data;

temp =table;
temptemp=table;
for isub=1:height(sdata)
    bdata = beh_data(beh_data.Subject==sdata.Subject(isub),:);
    temptemp = (bdata.Role==1).*(...
                                 (sdata.PTMzBuy(isub)>=0 ).*( (1-sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))+sdata.PTMzBuy(isub) ) + ...
                                 (sdata.PTMzBuy(isub)< 0 ).*  (1+sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))...
                                 ) + ...
               (bdata.Role==2).*(...
                                 (sdata.PTMzSell(isub)>=0).*( (1-sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))+sdata.PTMzSell(isub)) + ...
                                 (sdata.PTMzSell(isub)< 0).*  (1+sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))...
                                 ) ;
    temp=[temp;array2table(temptemp)];
end

beh_data=beh_data;
beh_data.PTMProbPred=temp.temptemp;

beh_data.PTMEndingPred=beh_data.PTMProbPred>=0.5;
beh_data.PTMChoicePred = (2-beh_data.Role) .* beh_data.PTMEndingPred + (beh_data.Role-1) .* (1-beh_data.PTMEndingPred);

for isub=1:height(sub_data)
    bdata=beh_data(beh_data.Subject==isub,:);
    sub_data.PTM3ChoicePred(isub)=mean(bdata.Choice==bdata.PTMChoicePred); 
    ChoicePred(isub,idx)=mean(bdata.Choice==bdata.PTMChoicePred); 
end   

%% PTM4: no role effect on intercept

% params(1): lottery weighting for buyer
% params(2): decisiveness
% params(3): response bias for buyer
% params(4): lottery weighting for seller
% params(5): response bias for seller
% params(6): fixed utility for buyer
% params(6): fixed utility for seller

idx=4;
for isub = 1:height(sub_data)

    theta = [1    1    0    1    0    0     ]; % Initial values
    lb =    [0    0   -1    0   -1  -Inf    ];
    ub =    [Inf Inf   1   Inf   1   Inf    ];    
    options = optimset('Display','iter','TolFun',1e-8,'MaxFunEvals',10000,'Display','notify','MaxIter',10000);   
    wdata = beh_data(beh_data.Subject==sub_data.Subject(isub),:);
    [MLE{idx}(isub,1:length(theta)),NegLogL(isub,idx),exitflag(isub,idx)] = fmincon(@ptm4,theta,[],[],[],[],lb,ub,[],options,wdata);
    LogL(isub,idx) = -NegLogL(isub,idx);
    [AIC(isub,idx),BIC(isub,idx)] = aicbic(LogL(isub,idx),length(theta),height(wdata));
    
end

sub_data.PTMwLotteryBuy=MLE{idx}(:,1);
sub_data.PTMwLotterySell=MLE{idx}(:,4);

sub_data.PTMzBuy=MLE{idx}(:,3);
sub_data.PTMzSell=MLE{idx}(:,5);

sub_data.PTMInterBuy=MLE{idx}(:,6);
sub_data.PTMInterSell=MLE{idx}(:,6);

sub_data.PTMMiuBuy=MLE{idx}(:,2);
sub_data.PTMMiuSell=MLE{idx}(:,2);

sub_data.PTMResBias=sub_data.PTMzSell-sub_data.PTMzBuy;
sub_data.PTMValBias=sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy;
sub_data.PTMInterBias=sub_data.PTMInterSell-sub_data.PTMInterBuy;
sub_data.PTMMiuBias=sub_data.PTMMiuSell./sub_data.PTMMiuBuy;

sub_data.PTMValBiasLog=log(sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy);

% table
A=mean([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
B=std([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
rowNames = {'mean','std'};
colNames = {'PTMwLotteryBuy','PTMwLotterySell','PTMValBias','PTMzBuy','PTMzSell','PTMResBias','PTMInterBuy','PTMInterSell','PTMInterBias','PTMMiuBuy','PTMMiuSell','PTMMiuBias','PTMValBiasLog'};
disp(array2table([A;B],'RowNames',rowNames,'VariableNames',colNames));
[A;B];

%% PTM1 predict

sdata=sub_data;

temp =table;
temptemp=table;
for isub=1:height(sdata)
    bdata = beh_data(beh_data.Subject==sdata.Subject(isub),:);
    temptemp = (bdata.Role==1).*(...
                                 (sdata.PTMzBuy(isub)>=0 ).*( (1-sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))+sdata.PTMzBuy(isub) ) + ...
                                 (sdata.PTMzBuy(isub)< 0 ).*  (1+sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))...
                                 ) + ...
               (bdata.Role==2).*(...
                                 (sdata.PTMzSell(isub)>=0).*( (1-sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))+sdata.PTMzSell(isub)) + ...
                                 (sdata.PTMzSell(isub)< 0).*  (1+sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))...
                                 ) ;
    temp=[temp;array2table(temptemp)];
end

beh_data=beh_data;
beh_data.PTMProbPred=temp.temptemp;

beh_data.PTMEndingPred=beh_data.PTMProbPred>=0.5;
beh_data.PTMChoicePred = (2-beh_data.Role) .* beh_data.PTMEndingPred + (beh_data.Role-1) .* (1-beh_data.PTMEndingPred);

for isub=1:height(sub_data)
    bdata=beh_data(beh_data.Subject==isub,:);
    sub_data.PTM4ChoicePred(isub)=mean(bdata.Choice==bdata.PTMChoicePred); 
    ChoicePred(isub,idx)=mean(bdata.Choice==bdata.PTMChoicePred); 
end   

%% PTM5: weighting=1

% 1: lottery weighting for buyer
% params(1): decisiveness
% params(2): response bias for buyer
% 1: lottery weighting for seller
% params(3): response bias for seller
% params(4): fixed utility for buyer
% params(5): fixed utility for seller

idx=5;
for isub = 1:height(sub_data)

    theta = [     1    0         0    0    0 ]; % Initial values
    lb =    [     0   -1        -1  -Inf -Inf];
    ub =    [    Inf   1         1   Inf  Inf];    
    options = optimset('Display','iter','TolFun',1e-8,'MaxFunEvals',10000,'Display','notify','MaxIter',10000);   
    wdata = beh_data(beh_data.Subject==sub_data.Subject(isub),:);
    [MLE{idx}(isub,1:length(theta)),NegLogL(isub,idx),exitflag(isub,idx)] = fmincon(@ptm5,theta,[],[],[],[],lb,ub,[],options,wdata);
    LogL(isub,idx) = -NegLogL(isub,idx);
    [AIC(isub,idx),BIC(isub,idx)] = aicbic(LogL(isub,idx),length(theta),height(wdata));
    
end

sub_data.PTMwLotteryBuy(:)=1;
sub_data.PTMwLotterySell(:)=1;

sub_data.PTMzBuy=MLE{idx}(:,2);
sub_data.PTMzSell=MLE{idx}(:,3);

sub_data.PTMInterBuy=MLE{idx}(:,4);
sub_data.PTMInterSell=MLE{idx}(:,5);

sub_data.PTMMiuBuy=MLE{idx}(:,1);
sub_data.PTMMiuSell=MLE{idx}(:,1);

sub_data.PTMResBias=sub_data.PTMzSell-sub_data.PTMzBuy;
sub_data.PTMValBias=sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy;
sub_data.PTMInterBias=sub_data.PTMInterSell-sub_data.PTMInterBuy;
sub_data.PTMMiuBias=sub_data.PTMMiuSell./sub_data.PTMMiuBuy;

sub_data.PTMValBiasLog=log(sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy);

% table
A=mean([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
B=std([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
rowNames = {'mean','std'};
colNames = {'PTMwLotteryBuy','PTMwLotterySell','PTMValBias','PTMzBuy','PTMzSell','PTMResBias','PTMInterBuy','PTMInterSell','PTMInterBias','PTMMiuBuy','PTMMiuSell','PTMMiuBias','PTMValBiasLog'};
disp(array2table([A;B],'RowNames',rowNames,'VariableNames',colNames));
[A;B];

%% PTM1 predict

sdata=sub_data;

temp =table;
temptemp=table;
for isub=1:height(sdata)
    bdata = beh_data(beh_data.Subject==sdata.Subject(isub),:);
    temptemp = (bdata.Role==1).*(...
                                 (sdata.PTMzBuy(isub)>=0 ).*( (1-sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))+sdata.PTMzBuy(isub) ) + ...
                                 (sdata.PTMzBuy(isub)< 0 ).*  (1+sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))...
                                 ) + ...
               (bdata.Role==2).*(...
                                 (sdata.PTMzSell(isub)>=0).*( (1-sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))+sdata.PTMzSell(isub)) + ...
                                 (sdata.PTMzSell(isub)< 0).*  (1+sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))...
                                 ) ;
    temp=[temp;array2table(temptemp)];
end

beh_data=beh_data;
beh_data.PTMProbPred=temp.temptemp;

beh_data.PTMEndingPred=beh_data.PTMProbPred>=0.5;
beh_data.PTMChoicePred = (2-beh_data.Role) .* beh_data.PTMEndingPred + (beh_data.Role-1) .* (1-beh_data.PTMEndingPred);

for isub=1:height(sub_data)
    bdata=beh_data(beh_data.Subject==isub,:);
    sub_data.PTM5ChoicePred(isub)=mean(bdata.Choice==bdata.PTMChoicePred); 
    ChoicePred(isub,idx)=mean(bdata.Choice==bdata.PTMChoicePred); 
end   

%% PTM6: response=0

% params(1): lottery weighting for buyer
% params(2): decisiveness
% 0: response bias for buyer
% params(3): lottery weighting for seller
% 0: response bias for seller
% params(4): fixed utility for buyer
% params(5): fixed utility for seller

idx=6;
for isub = 1:height(sub_data)

    theta = [1    1         1         0    0 ]; % Initial values
    lb =    [0    0         0       -Inf -Inf];
    ub =    [Inf Inf       Inf       Inf  Inf];    
    options = optimset('Display','iter','TolFun',1e-8,'MaxFunEvals',10000,'Display','notify','MaxIter',10000);   
    wdata = beh_data(beh_data.Subject==sub_data.Subject(isub),:);
    [MLE{idx}(isub,1:length(theta)),NegLogL(isub,idx),exitflag(isub,idx)] = fmincon(@ptm6,theta,[],[],[],[],lb,ub,[],options,wdata);
    LogL(isub,idx) = -NegLogL(isub,idx);
    [AIC(isub,idx),BIC(isub,idx)] = aicbic(LogL(isub,idx),length(theta),height(wdata));
    
end

sub_data.PTMwLotteryBuy=MLE{idx}(:,1);
sub_data.PTMwLotterySell=MLE{idx}(:,3);

sub_data.PTMzBuy(:)=0;
sub_data.PTMzSell(:)=0;

sub_data.PTMInterBuy=MLE{idx}(:,4);
sub_data.PTMInterSell=MLE{idx}(:,5);

sub_data.PTMMiuBuy=MLE{idx}(:,2);
sub_data.PTMMiuSell=MLE{idx}(:,2);

sub_data.PTMResBias=sub_data.PTMzSell-sub_data.PTMzBuy;
sub_data.PTMValBias=sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy;
sub_data.PTMInterBias=sub_data.PTMInterSell-sub_data.PTMInterBuy;
sub_data.PTMMiuBias=sub_data.PTMMiuSell./sub_data.PTMMiuBuy;

sub_data.PTMValBiasLog=log(sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy);

% table
A=mean([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
B=std([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
rowNames = {'mean','std'};
colNames = {'PTMwLotteryBuy','PTMwLotterySell','PTMValBias','PTMzBuy','PTMzSell','PTMResBias','PTMInterBuy','PTMInterSell','PTMInterBias','PTMMiuBuy','PTMMiuSell','PTMMiuBias','PTMValBiasLog'};
disp(array2table([A;B],'RowNames',rowNames,'VariableNames',colNames));
[A;B];

%% PTM1 predict

sdata=sub_data;

temp =table;
temptemp=table;
for isub=1:height(sdata)
    bdata = beh_data(beh_data.Subject==sdata.Subject(isub),:);
    temptemp = (bdata.Role==1).*(...
                                 (sdata.PTMzBuy(isub)>=0 ).*( (1-sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))+sdata.PTMzBuy(isub) ) + ...
                                 (sdata.PTMzBuy(isub)< 0 ).*  (1+sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))...
                                 ) + ...
               (bdata.Role==2).*(...
                                 (sdata.PTMzSell(isub)>=0).*( (1-sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))+sdata.PTMzSell(isub)) + ...
                                 (sdata.PTMzSell(isub)< 0).*  (1+sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))...
                                 ) ;
    temp=[temp;array2table(temptemp)];
end

beh_data=beh_data;
beh_data.PTMProbPred=temp.temptemp;

beh_data.PTMEndingPred=beh_data.PTMProbPred>=0.5;
beh_data.PTMChoicePred = (2-beh_data.Role) .* beh_data.PTMEndingPred + (beh_data.Role-1) .* (1-beh_data.PTMEndingPred);

for isub=1:height(sub_data)
    bdata=beh_data(beh_data.Subject==isub,:);
    sub_data.PTM6ChoicePred(isub)=mean(bdata.Choice==bdata.PTMChoicePred); 
    ChoicePred(isub,idx)=mean(bdata.Choice==bdata.PTMChoicePred); 
end   

%% PTM7: intercept=0

% params(1): lottery weighting for buyer
% params(2): decisiveness
% params(3): response bias for buyer
% params(4): lottery weighting for seller
% params(5): response bias for seller
% 0: fixed utility for buyer
% 0: fixed utility for seller

idx=7;
for isub = 1:height(sub_data)

    theta = [1    1    0    1    0           ]; % Initial values
    lb =    [0    0   -1    0   -1           ];
    ub =    [Inf Inf   1   Inf   1           ];    
    options = optimset('Display','iter','TolFun',1e-8,'MaxFunEvals',10000,'Display','notify','MaxIter',10000);   
    wdata = beh_data(beh_data.Subject==sub_data.Subject(isub),:);
    [MLE{idx}(isub,1:length(theta)),NegLogL(isub,idx),exitflag(isub,idx)] = fmincon(@ptm7,theta,[],[],[],[],lb,ub,[],options,wdata);
    LogL(isub,idx) = -NegLogL(isub,idx);
    [AIC(isub,idx),BIC(isub,idx)] = aicbic(LogL(isub,idx),length(theta),height(wdata));
    
end

sub_data.PTMwLotteryBuy=MLE{idx}(:,1);
sub_data.PTMwLotterySell=MLE{idx}(:,4);

sub_data.PTMzBuy=MLE{idx}(:,3);
sub_data.PTMzSell=MLE{idx}(:,5);

sub_data.PTMInterBuy(:)=0;
sub_data.PTMInterSell(:)=0;

sub_data.PTMMiuBuy=MLE{idx}(:,2);
sub_data.PTMMiuSell=MLE{idx}(:,2);

sub_data.PTMResBias=sub_data.PTMzSell-sub_data.PTMzBuy;
sub_data.PTMValBias=sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy;
sub_data.PTMInterBias=sub_data.PTMInterSell-sub_data.PTMInterBuy;
sub_data.PTMMiuBias=sub_data.PTMMiuSell./sub_data.PTMMiuBuy;

sub_data.PTMValBiasLog=log(sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy);

% table
A=mean([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
B=std([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
rowNames = {'mean','std'};
colNames = {'PTMwLotteryBuy','PTMwLotterySell','PTMValBias','PTMzBuy','PTMzSell','PTMResBias','PTMInterBuy','PTMInterSell','PTMInterBias','PTMMiuBuy','PTMMiuSell','PTMMiuBias','PTMValBiasLog'};
disp(array2table([A;B],'RowNames',rowNames,'VariableNames',colNames));
[A;B];

%% PTM1 predict

sdata=sub_data;

temp =table;
temptemp=table;
for isub=1:height(sdata)
    bdata = beh_data(beh_data.Subject==sdata.Subject(isub),:);
    temptemp = (bdata.Role==1).*(...
                                 (sdata.PTMzBuy(isub)>=0 ).*( (1-sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))+sdata.PTMzBuy(isub) ) + ...
                                 (sdata.PTMzBuy(isub)< 0 ).*  (1+sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))...
                                 ) + ...
               (bdata.Role==2).*(...
                                 (sdata.PTMzSell(isub)>=0).*( (1-sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))+sdata.PTMzSell(isub)) + ...
                                 (sdata.PTMzSell(isub)< 0).*  (1+sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))...
                                 ) ;
    temp=[temp;array2table(temptemp)];
end

beh_data=beh_data;
beh_data.PTMProbPred=temp.temptemp;

beh_data.PTMEndingPred=beh_data.PTMProbPred>=0.5;
beh_data.PTMChoicePred = (2-beh_data.Role) .* beh_data.PTMEndingPred + (beh_data.Role-1) .* (1-beh_data.PTMEndingPred);

for isub=1:height(sub_data)
    bdata=beh_data(beh_data.Subject==isub,:);
    sub_data.PTM7ChoicePred(isub)=mean(bdata.Choice==bdata.PTMChoicePred); 
    ChoicePred(isub,idx)=mean(bdata.Choice==bdata.PTMChoicePred); 
end   

%% PTM8: include role effect on decisiveness

% params(1): lottery weighting for buyer
% params(2): decisiveness
% params(3): response bias for buyer
% params(4): lottery weighting for seller
% params(5): response bias for seller
% params(6): fixed utility for buyer
% params(7): fixed utility for seller
% params(8): decisiveness for seller

idx=8;
for isub = 1:height(sub_data)

    theta = [1    1    0    1    0    0    0   1 ]; % Initial values
    lb =    [0    0   -1    0   -1  -Inf -Inf  0 ];
    ub =    [Inf Inf   1   Inf   1   Inf  Inf Inf];    
    options = optimset('Display','iter','TolFun',1e-8,'MaxFunEvals',10000,'Display','notify','MaxIter',10000);   
    wdata = beh_data(beh_data.Subject==sub_data.Subject(isub),:);
    [MLE{idx}(isub,1:length(theta)),NegLogL(isub,idx),exitflag(isub,idx)] = fmincon(@ptm8,theta,[],[],[],[],lb,ub,[],options,wdata);
    LogL(isub,idx) = -NegLogL(isub,idx);
    [AIC(isub,idx),BIC(isub,idx)] = aicbic(LogL(isub,idx),length(theta),height(wdata));
    
end

sub_data.PTMwLotteryBuy=MLE{idx}(:,1);
sub_data.PTMwLotterySell=MLE{idx}(:,4);

sub_data.PTMzBuy=MLE{idx}(:,3);
sub_data.PTMzSell=MLE{idx}(:,5);

sub_data.PTMInterBuy=MLE{idx}(:,6);
sub_data.PTMInterSell=MLE{idx}(:,7);

sub_data.PTMMiuBuy=MLE{idx}(:,2);
sub_data.PTMMiuSell=MLE{idx}(:,8);

sub_data.PTMResBias=sub_data.PTMzSell-sub_data.PTMzBuy;
sub_data.PTMValBias=sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy;
sub_data.PTMInterBias=sub_data.PTMInterSell-sub_data.PTMInterBuy;
sub_data.PTMMiuBias=sub_data.PTMMiuSell./sub_data.PTMMiuBuy;

sub_data.PTMValBiasLog=log(sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy);

% table
A=mean([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
B=std([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
rowNames = {'mean','std'};
colNames = {'PTMwLotteryBuy','PTMwLotterySell','PTMValBias','PTMzBuy','PTMzSell','PTMResBias','PTMInterBuy','PTMInterSell','PTMInterBias','PTMMiuBuy','PTMMiuSell','PTMMiuBias','PTMValBiasLog'};
disp(array2table([A;B],'RowNames',rowNames,'VariableNames',colNames));
[A;B];

%% PTM1 predict

sdata=sub_data;

temp =table;
temptemp=table;
for isub=1:height(sdata)
    bdata = beh_data(beh_data.Subject==sdata.Subject(isub),:);
    temptemp = (bdata.Role==1).*(...
                                 (sdata.PTMzBuy(isub)>=0 ).*( (1-sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))+sdata.PTMzBuy(isub) ) + ...
                                 (sdata.PTMzBuy(isub)< 0 ).*  (1+sdata.PTMzBuy(isub))./(1+exp(-sdata.PTMMiuBuy(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterBuy(isub)))...
                                 ) + ...
               (bdata.Role==2).*(...
                                 (sdata.PTMzSell(isub)>=0).*( (1-sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))+sdata.PTMzSell(isub)) + ...
                                 (sdata.PTMzSell(isub)< 0).*  (1+sdata.PTMzSell(isub))./(1+exp(-sdata.PTMMiuSell(isub)*(0.5.*sdata.PTMwLotteryBuy(isub).*bdata.Lottery-bdata.Price) + sdata.PTMInterSell(isub)))...
                                 ) ;
    temp=[temp;array2table(temptemp)];
end

beh_data=beh_data;
beh_data.PTMProbPred=temp.temptemp;

beh_data.PTMEndingPred=beh_data.PTMProbPred>=0.5;
beh_data.PTMChoicePred = (2-beh_data.Role) .* beh_data.PTMEndingPred + (beh_data.Role-1) .* (1-beh_data.PTMEndingPred);

for isub=1:height(sub_data)
    bdata=beh_data(beh_data.Subject==isub,:);
    sub_data.PTM8ChoicePred(isub)=mean(bdata.Choice==bdata.PTMChoicePred); 
    ChoicePred(isub,idx)=mean(bdata.Choice==bdata.PTMChoicePred); 
end   

%% PTM model comparison

sum(AIC)';
sum(BIC)';
sum(LogL)';
sum(ChoicePred)

%%









%% PTM1: re-save PTM estimated variables from PTM1

% params(1): lottery weighting for buyer
% params(2): decisiveness
% params(3): response bias for buyer
% params(4): lottery weighting for seller
% params(5): response bias for seller
% params(6): fixed utility for buyer
% params(7): fixed utility for seller

idx=1;
for isub = 1:height(sub_data)

    theta = [1    1    0    1    0    0    0 ]; % Initial values
    lb =    [0    0   -1    0   -1  -Inf -Inf];
    ub =    [Inf Inf   1   Inf   1   Inf  Inf];    
    options = optimset('Display','iter','TolFun',1e-8,'MaxFunEvals',10000,'Display','notify','MaxIter',10000);   
    wdata = beh_data(beh_data.Subject==sub_data.Subject(isub),:);
    [MLE{idx}(isub,1:length(theta)),NegLogL(isub,idx),exitflag(isub,idx)] = fmincon(@ptm1,theta,[],[],[],[],lb,ub,[],options,wdata);
    LogL(isub,idx) = -NegLogL(isub,idx);
    [AIC(isub,idx),BIC(isub,idx)] = aicbic(LogL(isub,idx),length(theta),height(wdata));
    
end

sub_data.PTMwLotteryBuy=MLE{idx}(:,1);
sub_data.PTMwLotterySell=MLE{idx}(:,4);

sub_data.PTMzBuy=MLE{idx}(:,3);
sub_data.PTMzSell=MLE{idx}(:,5);

sub_data.PTMInterBuy=MLE{idx}(:,6);
sub_data.PTMInterSell=MLE{idx}(:,7);

sub_data.PTMMiuBuy=MLE{idx}(:,2);
sub_data.PTMMiuSell=MLE{idx}(:,2);

sub_data.PTMResBias=sub_data.PTMzSell-sub_data.PTMzBuy;
sub_data.PTMValBias=sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy;
sub_data.PTMInterBias=sub_data.PTMInterSell-sub_data.PTMInterBuy;
sub_data.PTMMiuBias=sub_data.PTMMiuSell./sub_data.PTMMiuBuy;

sub_data.PTMValBiasLog=log(sub_data.PTMwLotterySell./sub_data.PTMwLotteryBuy);

% table
A=mean([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
B=std([sub_data.PTMwLotteryBuy,sub_data.PTMwLotterySell,sub_data.PTMValBias,sub_data.PTMzBuy,sub_data.PTMzSell,sub_data.PTMResBias,sub_data.PTMInterBuy,sub_data.PTMInterSell,sub_data.PTMInterBias,sub_data.PTMMiuBuy,sub_data.PTMMiuSell,sub_data.PTMMiuBias,sub_data.PTMValBiasLog]);
rowNames = {'mean','std'};
colNames = {'PTMwLotteryBuy','PTMwLotterySell','PTMValBias','PTMzBuy','PTMzSell','PTMResBias','PTMInterBuy','PTMInterSell','PTMInterBias','PTMMiuBuy','PTMMiuSell','PTMMiuBias','PTMValBiasLog'};
disp(array2table([A;B],'RowNames',rowNames,'VariableNames',colNames));
[A;B];

%% 
writetable(sub_data,[root_path,'\data\subject\subject_ptm.csv']);
toc;
